Access to care

library(tidyverse)
library(rlang)
library(leaflet)

hospitals <- read_csv("../data/hospital_list.csv")

Parsed with column specification:
cols(
  .default = col_character(),
  zip = col_double(),
  `Phone Number` = col_double(),
  score = col_double(),
  longitude = col_double(),
  latitude = col_double(),
  population = col_double(),
  poverty = col_double(),
  percent = col_double(),
  poverty_percent = col_double()
)
See spec(...) for full column specifications.

hospitals

hospital_locations <- hospitals  %>%
  select(state, longitude, latitude)
write_rds(hospital_locations, "../data/hospital_locations.rds")
hospital_locations

population <- usmap::countypop

hospital_count <- hospitals %>%
  count(state, original_county) 


state_hospitals <- population %>%
  left_join(hospital_count, by = c("abbr" = "state", "county" = "original_county")) %>%
  mutate(hospitals = ifelse(is.na(n), 0, n)) %>%
  select(
    fips, 
    state = abbr, 
    county, 
    population = pop_2015,
    hospitals
    )

state_hospitals

set.seed(100)

hospital_sample <- state_hospitals %>%
  sample_frac(0.3)

hospital_sample

model <- lm(hospitals ~ population, data = hospital_sample)
write_rds(model, "../data/model.rds")
summary(model)

Call:
lm(formula = hospitals ~ population, data = hospital_sample)

Residuals:
     Min       1Q   Median       3Q      Max 
-10.7909  -0.7106   0.0908   0.2498   8.3805 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 7.027e-01  3.942e-02   17.83   <2e-16 ***
population  7.907e-06  1.182e-07   66.90   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.149 on 941 degrees of freedom
Multiple R-squared:  0.8263,    Adjusted R-squared:  0.8261 
F-statistic:  4475 on 1 and 941 DF,  p-value: < 2.2e-16

predictions <- predict(model, 
  newdata = state_hospitals, 
  interval = "prediction") %>%
  as_tibble() %>%
  mutate_all(round)

predictions

hospital_results <- state_hospitals %>%
  bind_cols(predictions) %>%
  mutate(result = case_when(
    hospitals < lwr ~ -1,
    hospitals > upr ~ 1,
    TRUE ~ 0)) %>%
  mutate(
    county_key = str_remove_all(county, " County"),
    county_key = str_remove_all(county_key, " Parish"),
    county_key = str_remove_all(county_key, " city"),
    county_key = str_to_lower(county_key),
    county_key = str_remove_all(county_key, " "),
    county_key = str_replace_all(county_key, "st. ", "saint"),
      )  

write_rds(hospital_results, "../data/hospitals.rds")

hospital_results 

shapes <- read_rds("../data/shapes.rds") %>%
  mutate(
    county_key = str_to_lower(county),
    county_key = str_remove_all(county_key, " "),
    county_key = str_replace_all(county_key, "st. ", "saint")
  ) 

shapes 

county_hospitals <- hospital_results  %>%
  select(- county) %>%
  left_join(shapes, by = c("state", "county_key")) %>%
  filter(state != "PR", state != "AK") %>%
  filter(!is.na(county))

county_hospitals %>%
  mutate(popup = paste0("County:", county, " Population:", population, " Hospitals:", hospitals)) %>%
  select(popup, everything())

states <- county_hospitals %>%
  group_by(state) %>%
  summarise() %>%
  pull()

under <- "#CC79A7"
over <- "#0072B2"
at_level <- "#008b00"
hospital_color <- "#F0E442"
st <- c("LA")
#st <- states
include <- c(-1, 1, 0)
add_hospitals <- 0

initial_map <- leaflet() %>%
  addProviderTiles(providers$CartoDB) %>%
  addLegend("bottomright", 
            color  = c(under,over ,at_level, hospital_color), 
            labels = c("Less hopitals than expected",
                       "More hospitals than expected", "Within Range", 
                       "Hospital Location"), 
            title  = "Legend",opacity = 0.5)

locations <- hospital_locations %>%
  filter(state %in% st) %>%
  select(longitude, latitude) %>%
  mutate_all(round, 3) %>%
  count(longitude, latitude) %>%
  mutate(popup = paste0("Hospitals: ", n))

counties <- county_hospitals %>% 
  filter(state %in% st,
         result %in% include) %>%
  mutate(popup = paste0("<b>", county, "</b>",
                        "<br/>Population: ", prettyNum(population, big.mark = ","), 
                        "<br/>Hospitals: ", hospitals,
                        "<br/>Expected: ", fit
                        )) %>%  
  mutate(color = case_when(
    result ==  0  ~ at_level,
    result ==  1  ~ over,
    result == -1  ~ under
  )) 

county_map <- counties %>%
  transpose() %>%
  map(~{
    county <- .x
     transpose(county$data)  %>%
       map(~list(
         data = .x$data, 
         color = county$color,
         popup = county$popup
         ))
    }) %>%
  flatten() %>%
  reduce(
    ~ addPolygons(.x, 
                  lng = .y$data$long, 
                  lat = .y$data$lat, 
                  color = .y$color, 
                  popup = .y$popup,
                  weight = 1, fillOpacity = 0.3), 
    .init = initial_map)  


if(add_hospitals == 1) {
  county_map <- county_map %>%
    addCircleMarkers(
      lng = locations$longitude, 
      lat = locations$latitude,
      radius = 1.5 *locations$n,
      popup = locations$popup,
      fillColor = hospital_color, color = "gray", 
      fillOpacity = 0.7,weight = 1)
  } 

county_map

NA

LS0tDQp0aXRsZTogIkFjY2VzcyB0byBjYXJlIg0Kb3V0cHV0OiBodG1sX25vdGVib29rDQotLS0NCg0KYGBge3J9DQpsaWJyYXJ5KHRpZHl2ZXJzZSkNCmxpYnJhcnkocmxhbmcpDQpsaWJyYXJ5KGxlYWZsZXQpDQpgYGANCg0KYGBge3J9DQpob3NwaXRhbHMgPC0gcmVhZF9jc3YoIi4uL2RhdGEvaG9zcGl0YWxfbGlzdC5jc3YiKQ0KaG9zcGl0YWxzDQpgYGANCg0KYGBge3J9DQpob3NwaXRhbF9sb2NhdGlvbnMgPC0gaG9zcGl0YWxzICAlPiUNCiAgc2VsZWN0KHN0YXRlLCBsb25naXR1ZGUsIGxhdGl0dWRlKQ0Kd3JpdGVfcmRzKGhvc3BpdGFsX2xvY2F0aW9ucywgIi4uL2RhdGEvaG9zcGl0YWxfbG9jYXRpb25zLnJkcyIpDQpob3NwaXRhbF9sb2NhdGlvbnMNCmBgYA0KDQpgYGB7cn0NCnBvcHVsYXRpb24gPC0gdXNtYXA6OmNvdW50eXBvcA0KDQpob3NwaXRhbF9jb3VudCA8LSBob3NwaXRhbHMgJT4lDQogIGNvdW50KHN0YXRlLCBvcmlnaW5hbF9jb3VudHkpIA0KDQoNCnN0YXRlX2hvc3BpdGFscyA8LSBwb3B1bGF0aW9uICU+JQ0KICBsZWZ0X2pvaW4oaG9zcGl0YWxfY291bnQsIGJ5ID0gYygiYWJiciIgPSAic3RhdGUiLCAiY291bnR5IiA9ICJvcmlnaW5hbF9jb3VudHkiKSkgJT4lDQogIG11dGF0ZShob3NwaXRhbHMgPSBpZmVsc2UoaXMubmEobiksIDAsIG4pKSAlPiUNCiAgc2VsZWN0KA0KICAgIGZpcHMsIA0KICAgIHN0YXRlID0gYWJiciwgDQogICAgY291bnR5LCANCiAgICBwb3B1bGF0aW9uID0gcG9wXzIwMTUsDQogICAgaG9zcGl0YWxzDQogICAgKQ0KDQpzdGF0ZV9ob3NwaXRhbHMNCmBgYA0KDQoNCg0KDQpgYGB7cn0NCnNldC5zZWVkKDEwMCkNCg0KaG9zcGl0YWxfc2FtcGxlIDwtIHN0YXRlX2hvc3BpdGFscyAlPiUNCiAgc2FtcGxlX2ZyYWMoMC4zKQ0KDQpob3NwaXRhbF9zYW1wbGUNCmBgYA0KDQpgYGB7cn0NCm1vZGVsIDwtIGxtKGhvc3BpdGFscyB+IHBvcHVsYXRpb24sIGRhdGEgPSBob3NwaXRhbF9zYW1wbGUpDQp3cml0ZV9yZHMobW9kZWwsICIuLi9kYXRhL21vZGVsLnJkcyIpDQpzdW1tYXJ5KG1vZGVsKQ0KYGBgDQoNCmBgYHtyfQ0KcHJlZGljdGlvbnMgPC0gcHJlZGljdChtb2RlbCwgDQogIG5ld2RhdGEgPSBzdGF0ZV9ob3NwaXRhbHMsIA0KICBpbnRlcnZhbCA9ICJwcmVkaWN0aW9uIikgJT4lDQogIGFzX3RpYmJsZSgpICU+JQ0KICBtdXRhdGVfYWxsKHJvdW5kKQ0KDQpwcmVkaWN0aW9ucw0KYGBgDQoNCmBgYHtyfQ0KaG9zcGl0YWxfcmVzdWx0cyA8LSBzdGF0ZV9ob3NwaXRhbHMgJT4lDQogIGJpbmRfY29scyhwcmVkaWN0aW9ucykgJT4lDQogIG11dGF0ZShyZXN1bHQgPSBjYXNlX3doZW4oDQogICAgaG9zcGl0YWxzIDwgbHdyIH4gLTEsDQogICAgaG9zcGl0YWxzID4gdXByIH4gMSwNCiAgICBUUlVFIH4gMCkpICU+JQ0KICBtdXRhdGUoDQogICAgY291bnR5X2tleSA9IHN0cl9yZW1vdmVfYWxsKGNvdW50eSwgIiBDb3VudHkiKSwNCiAgICBjb3VudHlfa2V5ID0gc3RyX3JlbW92ZV9hbGwoY291bnR5X2tleSwgIiBQYXJpc2giKSwNCiAgICBjb3VudHlfa2V5ID0gc3RyX3JlbW92ZV9hbGwoY291bnR5X2tleSwgIiBjaXR5IiksDQogICAgY291bnR5X2tleSA9IHN0cl90b19sb3dlcihjb3VudHlfa2V5KSwNCiAgICBjb3VudHlfa2V5ID0gc3RyX3JlbW92ZV9hbGwoY291bnR5X2tleSwgIiAiKSwNCiAgICBjb3VudHlfa2V5ID0gc3RyX3JlcGxhY2VfYWxsKGNvdW50eV9rZXksICJzdC4gIiwgInNhaW50IiksDQogICAgICApICANCg0Kd3JpdGVfcmRzKGhvc3BpdGFsX3Jlc3VsdHMsICIuLi9kYXRhL2hvc3BpdGFscy5yZHMiKQ0KDQpob3NwaXRhbF9yZXN1bHRzIA0KYGBgDQoNCiMjIFN0YXRlIG1hcA0KDQpgYGB7cn0NCnNoYXBlcyA8LSByZWFkX3JkcygiLi4vZGF0YS9zaGFwZXMucmRzIikgJT4lDQogIG11dGF0ZSgNCiAgICBjb3VudHlfa2V5ID0gc3RyX3RvX2xvd2VyKGNvdW50eSksDQogICAgY291bnR5X2tleSA9IHN0cl9yZW1vdmVfYWxsKGNvdW50eV9rZXksICIgIiksDQogICAgY291bnR5X2tleSA9IHN0cl9yZXBsYWNlX2FsbChjb3VudHlfa2V5LCAic3QuICIsICJzYWludCIpDQogICkgDQoNCnNoYXBlcyANCmBgYA0KYGBge3J9DQpjb3VudHlfaG9zcGl0YWxzIDwtIGhvc3BpdGFsX3Jlc3VsdHMgICU+JQ0KICBzZWxlY3QoLSBjb3VudHkpICU+JQ0KICBsZWZ0X2pvaW4oc2hhcGVzLCBieSA9IGMoInN0YXRlIiwgImNvdW50eV9rZXkiKSkgJT4lDQogIGZpbHRlcihzdGF0ZSAhPSAiUFIiLCBzdGF0ZSAhPSAiQUsiKSAlPiUNCiAgZmlsdGVyKCFpcy5uYShjb3VudHkpKQ0KDQpjb3VudHlfaG9zcGl0YWxzICU+JQ0KICBtdXRhdGUocG9wdXAgPSBwYXN0ZTAoIkNvdW50eToiLCBjb3VudHksICIgUG9wdWxhdGlvbjoiLCBwb3B1bGF0aW9uLCAiIEhvc3BpdGFsczoiLCBob3NwaXRhbHMpKSAlPiUNCiAgc2VsZWN0KHBvcHVwLCBldmVyeXRoaW5nKCkpDQpgYGANCg0KYGBge3J9DQpzdGF0ZXMgPC0gY291bnR5X2hvc3BpdGFscyAlPiUNCiAgZ3JvdXBfYnkoc3RhdGUpICU+JQ0KICBzdW1tYXJpc2UoKSAlPiUNCiAgcHVsbCgpDQoNCmBgYA0KDQoNCmBgYHtyfQ0KDQp1bmRlciA8LSAiI0NDNzlBNyINCm92ZXIgPC0gIiMwMDcyQjIiDQphdF9sZXZlbCA8LSAiIzAwOGIwMCINCmhvc3BpdGFsX2NvbG9yIDwtICIjRjBFNDQyIg0Kc3QgPC0gYygiTEEiKQ0KI3N0IDwtIHN0YXRlcw0KaW5jbHVkZSA8LSBjKC0xLCAxLCAwKQ0KYWRkX2hvc3BpdGFscyA8LSAwDQoNCmluaXRpYWxfbWFwIDwtIGxlYWZsZXQoKSAlPiUNCiAgYWRkUHJvdmlkZXJUaWxlcyhwcm92aWRlcnMkQ2FydG9EQikgJT4lDQogIGFkZExlZ2VuZCgiYm90dG9tcmlnaHQiLCANCiAgICAgICAgICAgIGNvbG9yICA9IGModW5kZXIsb3ZlciAsYXRfbGV2ZWwsIGhvc3BpdGFsX2NvbG9yKSwgDQogICAgICAgICAgICBsYWJlbHMgPSBjKCJMZXNzIGhvcGl0YWxzIHRoYW4gZXhwZWN0ZWQiLA0KICAgICAgICAgICAgICAgICAgICAgICAiTW9yZSBob3NwaXRhbHMgdGhhbiBleHBlY3RlZCIsICJXaXRoaW4gUmFuZ2UiLCANCiAgICAgICAgICAgICAgICAgICAgICAgIkhvc3BpdGFsIExvY2F0aW9uIiksIA0KICAgICAgICAgICAgdGl0bGUgID0gIkxlZ2VuZCIsb3BhY2l0eSA9IDAuNSkNCg0KbG9jYXRpb25zIDwtIGhvc3BpdGFsX2xvY2F0aW9ucyAlPiUNCiAgZmlsdGVyKHN0YXRlICVpbiUgc3QpICU+JQ0KICBzZWxlY3QobG9uZ2l0dWRlLCBsYXRpdHVkZSkgJT4lDQogIG11dGF0ZV9hbGwocm91bmQsIDMpICU+JQ0KICBjb3VudChsb25naXR1ZGUsIGxhdGl0dWRlKSAlPiUNCiAgbXV0YXRlKHBvcHVwID0gcGFzdGUwKCJIb3NwaXRhbHM6ICIsIG4pKQ0KDQpjb3VudGllcyA8LSBjb3VudHlfaG9zcGl0YWxzICU+JSANCiAgZmlsdGVyKHN0YXRlICVpbiUgc3QsDQogICAgICAgICByZXN1bHQgJWluJSBpbmNsdWRlKSAlPiUNCiAgbXV0YXRlKHBvcHVwID0gcGFzdGUwKCI8Yj4iLCBjb3VudHksICI8L2I+IiwNCiAgICAgICAgICAgICAgICAgICAgICAgICI8YnIvPlBvcHVsYXRpb246ICIsIHByZXR0eU51bShwb3B1bGF0aW9uLCBiaWcubWFyayA9ICIsIiksIA0KICAgICAgICAgICAgICAgICAgICAgICAgIjxici8+SG9zcGl0YWxzOiAiLCBob3NwaXRhbHMsDQogICAgICAgICAgICAgICAgICAgICAgICAiPGJyLz5FeHBlY3RlZDogIiwgZml0DQogICAgICAgICAgICAgICAgICAgICAgICApKSAlPiUgIA0KICBtdXRhdGUoY29sb3IgPSBjYXNlX3doZW4oDQogICAgcmVzdWx0ID09ICAwICB+IGF0X2xldmVsLA0KICAgIHJlc3VsdCA9PSAgMSAgfiBvdmVyLA0KICAgIHJlc3VsdCA9PSAtMSAgfiB1bmRlcg0KICApKSANCg0KY291bnR5X21hcCA8LSBjb3VudGllcyAlPiUNCiAgdHJhbnNwb3NlKCkgJT4lDQogIG1hcCh+ew0KICAgIGNvdW50eSA8LSAueA0KICAgICB0cmFuc3Bvc2UoY291bnR5JGRhdGEpICAlPiUNCiAgICAgICBtYXAofmxpc3QoDQogICAgICAgICBkYXRhID0gLngkZGF0YSwgDQogICAgICAgICBjb2xvciA9IGNvdW50eSRjb2xvciwNCiAgICAgICAgIHBvcHVwID0gY291bnR5JHBvcHVwDQogICAgICAgICApKQ0KICAgIH0pICU+JQ0KICBmbGF0dGVuKCkgJT4lDQogIHJlZHVjZSgNCiAgICB+IGFkZFBvbHlnb25zKC54LCANCiAgICAgICAgICAgICAgICAgIGxuZyA9IC55JGRhdGEkbG9uZywgDQogICAgICAgICAgICAgICAgICBsYXQgPSAueSRkYXRhJGxhdCwgDQogICAgICAgICAgICAgICAgICBjb2xvciA9IC55JGNvbG9yLCANCiAgICAgICAgICAgICAgICAgIHBvcHVwID0gLnkkcG9wdXAsDQogICAgICAgICAgICAgICAgICB3ZWlnaHQgPSAxLCBmaWxsT3BhY2l0eSA9IDAuMyksIA0KICAgIC5pbml0ID0gaW5pdGlhbF9tYXApICANCg0KDQppZihhZGRfaG9zcGl0YWxzID09IDEpIHsNCiAgY291bnR5X21hcCA8LSBjb3VudHlfbWFwICU+JQ0KICAgIGFkZENpcmNsZU1hcmtlcnMoDQogICAgICBsbmcgPSBsb2NhdGlvbnMkbG9uZ2l0dWRlLCANCiAgICAgIGxhdCA9IGxvY2F0aW9ucyRsYXRpdHVkZSwNCiAgICAgIHJhZGl1cyA9IDEuNSAqbG9jYXRpb25zJG4sDQogICAgICBwb3B1cCA9IGxvY2F0aW9ucyRwb3B1cCwNCiAgICAgIGZpbGxDb2xvciA9IGhvc3BpdGFsX2NvbG9yLCBjb2xvciA9ICJncmF5IiwgDQogICAgICBmaWxsT3BhY2l0eSA9IDAuNyx3ZWlnaHQgPSAxKQ0KICB9IA0KDQpjb3VudHlfbWFwDQogIA0KYGBgDQoNCg0KDQoNCg0KDQoNCg0KDQoNCg0KDQoNCg0KDQoNCg0KDQoNCg0KDQoNCg0KDQoNCg0KDQoNCg==